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We develop a microscopic magnetoelectric coupling in M3V2O8 (NVO) which gives rise to the 
trilinear phenomenological coupling used previously to explain the phase transition in which mag- 
netic and ferroelectric order parameters appear simultaneously. Using combined neutron scattering 
measurements and first-principles calculations of the phonons in NVO, we determine eleven phonons 
which can induce the observed spontaneous polarization. Among these eleven phonons, we find that 
a few of them can actually induce a significant dipole moment. Using the calculated atomic charges, 
we find that the required distortion to induce the observed dipole moment is very small ( 0.001 
A) and therefore it would be very difficult to observe the distortion by neutron-powder diffraction. 
Finally, we identify the derivatives of the exchange tensor with respect to atomic displacements 
which are needed for a microscopic model of a spin-phonon coupling in NVO and which we hope 
will be obtained from a fundamental quantum calculation such as LDA+U. We also analyze two 
toy models to illustrate that the Dzyaloskinskii-Moriya interaction is very important for coexisting 
of magnetic and ferroelectric order but it is not the only mechanism when the local site symmetry 
of the system is low enough. 

PACS numbers: 75.25.+Z, 75.10.Jm, 75.40.Gb 



I. INTRODUCTION 

Recently studies have identified a family of multifer- 
roics which display a phase transition in which there si- 
multaneously develops long-range incommensurate mag- 
netic and uniform ferroelectric order. Perhaps the most 
detailed studies have been carried out on the systems 
Ni 3 V 2 8 (NVmiS^ and TbMn0 3 (TMO)tf (For a re- 
view, see Ref. 0.) This phenomenon has been explained 3 
on the basis of a phenomenological model which invokes 
a Landau expansion in terms of the order parameters 
describing the incommensurate magnetic order and the 
order parameter describing the uniform spontaneous po- 
larization. Already from this treatment it was clear that 
a microscopic model would have to involve a trilinear 
interaction Hamiltonian proportional to the product of 
two spin variables and one displacement variable. Fur- 
thermore, the symmetry requirements of the phenomeno- 
logical model would naturally be realized by a proper mi- 
croscopic model. Accordingly, in this paper we present 
a detailed combined neutron scattering study and first- 
principles calculations of the optical phonons of NVO and 
thereby identify those having the right symmetry to in- 
duce a dipole moment. We also consider the expansion of 
the exchange tensor to first order in the generalized dis- 
placement coordinates. The aim of this programme is to 
determine which of the generalized displacements are rel- 
evant and which corresponding elements of the exchange 
tensor are needed for this microscopic calculation. This 
paper will therefore set the stage for a separate quantum 
calculation of the exchange tensor and its derivatives with 
respect to atomic displacements. On general grounds, 
one might expect the Dzyaloshinskii-Moriya 8 ' 9 (DM) in- 



teraction to play an important role in this calculation, 10 
and indeed, we will find this to be the case. However, 
the conclusion is a bit more subtle, in that displacement 
derivatives of many elements of the exchange tensor (and 
not just the antisymmetric one of the DM interaction) are 
needed for the microscopic calculation. We will show that 
simplistic models are somewhat misleading in that they 
may lead one to believe that the only exchange tensor ele- 
ments whose derivatives are important are the DM ones. 
Although our calculations may seem complicated, they 
are essential if one wishes to actually relate the magne- 
toelectric coupling to detailed quantum calculations (via 
LDA 11 or LDA+C/— schemes) of the strain-dependence 
of the exchange tensors. The methodology of the present 
paper can be extended in a straightforward way to TMO, 
for instance. 

Briefly this paper is organized as follows. In Sec. II 
we give an overview of the calculation. In Sec. Ill we 
discuss the first-principles calculations of the zone-center 
phonons and identify those phonons which transform like 
a vector and which are thus candidates to produce a spon- 
taneous polarization. In this section we also present the 
neutron scattering measurements of the phonon density 
of states (DOS), which is found to be in good agree- 
ment with the calculated spectrum. Next in Sec. IV we 
summarize the results of the determination of the spin 
structure^*^ In Sec. V we then use the symmetry opera- 
tions of the crystal to show how the phonon derivatives 
of the various exchange tensors in the unit cell are re- 
lated to one another. Then in Sec. VI we show that a 
mean-field treatment of this spin-phonon coupling leads 
to the results obtained previously^ in a phenomenolog- 
ical model. Here we give expressions for the spontaneous 
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polarization in terms of graidents of the excahnge ten- 
sor whose evaluation remains to be done. It would be 
nice to have a simple model to illustrate these results. 
However, our studies of two "toy models" in Sec. VII in- 
dicate that they do not reproduce some essential features 
of our complete calculation. Finally, our conclusions are 
summarized in Sec. VIII. 



II. OVERVIEW OF CALCULATION 

Here we give a brief qualitative overview of the calcu- 
lation. First we review the symmetry of the space group 
of NVO, Cmca (No. 64 in Ref. M). The space group 
operations (apart from primitive translations) are spec- 
ified in Table (Here and below, sites within the unit 
cell are given as fractions of the sides of the conventional 
unit cell, so that (x,y,z) denotes (xa,yb, zc).) We now 
describe the sets of crystallographically equivalent sites 
which the various atoms occupy. (Such a set of crystal- 
lographically equivalent sites is called a Wyckoff orbit.) 
There are six such orbits as shown in Table [n] The first 
two are those of the Ni atoms, the first consisting of the 
two Ni(l) (a) sites (which we call "cross-tie" sites) and 
the second consisting of the four Ni(2) (e) sites (which 
we call "spine" sites). The four V (f) sites comprise the 
third orbit and the oxygen sites are distributed into two 
(f) orbits, one containing four 0(1) atoms, the other con- 
taining four 0(2) atoms, and a (g) orbit containing eight 
0(3) atoms. (The letters a, e, f, g classi fy th e site symme- 
try according to the convention of Ref. Il3t The number 
of sites in the orbit as listed in Ref. [l3| is twice what 
we give here because here we consider the primitive unit 
cell rather than the conventional unit cell.) The loca- 
tions of these sites are specified in the second column of 
Table [H] Note that there are two formula units of NVO 
per unit cell. The Ni sites form buckled planes which 
resemble a kagome lattice and two such adjacent planes 
are shown in Fig. 1 where the buckling is omitted for 
simplicity. There one sees that the Ni spine sites form 
chains along the rr-direction between which the Ni cross- 
tie sites are situated, with bonds to nearest-neighboring 
spine sites which form a cross tie. To illustrate the use 
of Table [H] we find that the eight operations of Table |I] 
acting on (0,0,0) generate four copies of each of the 2 
a sites which are at (0,0,0) and (1/2,0,1/2). Similarly 
one can generate the eight g sites by applying in turn 
the eight operations of Table[I]to the site at (x, y, z). (In 
each case it may be necessary to bring the site back into 
the original unit cell via a primitive translation vector.) 

We now review briefly the nature of the ordered phases 
which occur as the temperature T is lowered at zero ex- 
ternal magnetic fieldiS^ At high temperatures the sys- 
tem is paramagnetic. When T is lowered through the 
value Tph ~ 9. IK, an incommensurate phase appears 
(which we call the high-temperature incommensurate or 
HTI phase) in which the Ni spins on the spine chains 
are oriented very nearly along the x-axis with a modula- 



Er = (x,y,z) 2 z r = (x + 1/2, y, z + 1/2) 

2 y r = (x + l/2,y,z + l/2) 2 x v = (x,y,z) 

Ir = (x,y,z) m z r = (x + \,y,z + |) 
m v r = (x + j, y, z + 1/2) m x r = (x,y,z) 

TABLE I: General positions within the primitive unit cell for 
Cmca which describe the symmetry operations of this space 
group. 2 a is a two-fold rotation (or screw) axis and m a is 
a mirror (or glide) plane. The primitive translation vectors 
are ai = (a/2)i + (b/2)j, a2 = (a/2)i — (b/2)j, and a3 = ck, 
where a = 5.92170A, 6 = 11.37105A, and c = 8.22638Ap^ 



Atoms 


(x/a,y/b,z/c) 


Wyckoff 


Decomposition 


Ni(l) 


(0,0,0) 


2a 


A u +2Bi«+2B2u+B3 M 


Ni(2) 


(l/4,y/b,l/4) 
y=0.1298 (0.1304) 


4e 


A„+2Bi u -|-B2ii-l-2B3„ 
A 9 +2Bi g +B 29 +2B 39 


V(l) 


(0,y,z) 
y=0.3762 (0.3762) 
z=0.1197 (0.1196) 


4f 


AU+2B1U+2B2U+B31, 
2A 9 +2Bi 9 +B2g-l-2B3 9 


0(1) 


(0,y/b,z/c) 
y=0.2481 (0.2490) 
z=0.2308 (0.2301) 


4f 


A u +2Bi u +2B2u-l-B3„ 
2A 9 +2Bi 9 +B2 9 +2B3 9 


0(2) 


(0,y/b,z/c) 
y=0.0011(0.0008) 
z=0.2444 (0.2441) 


4f 


A u +2Bi u +2B2 U +^3u 
2A 9 +2Bi 9 +B2 9 -l-2B3 9 


0(3) 


(x/a,y/b,z/c) 
x=0.2656(0.2703) 
y=0.1192 (0.1184) 
z=0.0002 (0.0012) 


8g 


3A„+3Bi u +3B2u+3B 3u 
3A 9 +3Bi 9 +3B 29 +3B 39 



TABLE II: In column 2 we give the Wyckoff positions, the 
fractional coordinates (in coumn 2) and their multiplicity (in 
column 3) of atoms in the Wyckoff orbits of NVO listed in 
column 1. In column 2 we give the values of the structural 
parameters (e.g. x, y, and z) as deduced from diffraction 
data-*— and the corresponding values we find from the struc- 
tural minimization are given in parentheses. In column 4 
we give the symmetry decomposition of the displacements of 
atoms in each of the Wyckoff orbits. 



tion vector also along i. (The axes are denoted either a, 
b, and c, or x, y, and z and corresponding unit vectors 
are denoted i, j, and k.) As the temperature is further 
lowered through the value Tjjl ~ 6.3K transverse or- 
der appears at the same incommensurate wavevector and 
also order appears on the cross-tie sites, as shown in Fig. 
1. We call this phase the low-temperature incommen- 
surate or LTI phase. Within experimental uncertainty, 
these two ordering transitions are continuous. As the 
temperature is lowered through the value Tlc ~ 4K, a 
discontinuous transition occurs in which a commensurate 
antiferromagnetic phase appears. In this phase antiferro- 
magnetism results from the arrangement of spins within 
the unit cell in such a way that the magnetic unit cell 
remains identical to the paramagnetic unit cell. 
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FIG. 1: Schematic diagram showing the x- and y-components 
of the spins in the LTI phase. We used the parameters: q = 
0.28(27r/a), a s , x = 1.6, a c , y = 1.4, b s , y = 1.3, b c , x = -2.2, 
and 4>l = <f>H + 7r/2. The small z-components of spin are 
not represented. The planes are buckled, so that alternately 
spine chains are displaced above and below the planes shown 
(but this buckling is not shown). In the HTI phase the cross- 
ties have negligible moments and the spine chains have the 
incommensurately modulated longitudinal moments similar 
to those shown. The a spin components are odd under 2 X 
and the b spin components are even under 2 X , where 2 X is a 
two-fold rotation about the x axis. 

Perhaps surprisingly, it was found that ferroelectricity 
coincides with the existence of LTI order, and this behav- 
ior was explained by a phcnomcnological model based on 
a Landau expansion in powers of order parameters de- 
scribing ferroelectricity and those needed to describe the 
magnetic ordering of the HTI and LTI phases £ The phe- 
nomenological interaction V between the magnetic and 
ferroelectric order parameters was given as 

V = a a,x,Ycr x (<l)(7Y(-q)P a , (1) 

a X,Y=L,H 

where <rx(q) is the order parameter describing incom- 
mensurate order at wavevector q characteristic of the 
HTI phase (for X=H) or for the additional ordering ap- 
pearing in the LTI phase (for X=L) and a labels the 
Cartesian component of the uniform spontaneous polar- 
ization vector P. Using the symmetry properties of the 
order parameters o\f/(q) and ox(q) it was shown that 
only o>b,L,H and cib,H,L are nonzero. This result provides a 
phenomenological explanation for the experimental find- 
ing that a nonzero polarization is induced by incommen- 
surate magnetism only in the LTI phase and then only 
with P along the b axis. 

This phenomenological theory elucidates the symme- 
try of the magnetically induced ferroelectric state. To 
develop an analogous microscopic theory, one could treat 
a Hubbard-like Hamiltonian involving the 3d electrons of 



the Ni ions and also the 2p electrons of the O ions. From 
such a treatment one can obtain the spin Hamiltonian 
which describes the "low energy" sector of this Hubbard- 
like Hamiltonian. This low-energy sector is obtained by 
removing states whose energy relative to the low energy 
sector involves Coulomb integrals. In this approach one 
develops a canonical transformation to eliminate hopping 
matrix elements. This approach is obviously most ap- 
propriate for insulators. This type of calculation was 
formalized by Anderson^ who thereby obtained the an- 
tiferromagnetic Heisenberg model from a superexchangc 
mechanism. More recently this approach has been ap- 
plied to LaCu 2 0^ii and LaTi0 3 i2iia Here we will use 
the phenomenological model to deduce the form of the 
low energy spin Hamiltonian which describes magneti- 
cally induced ferroelectricity in NVO. We are currently 
analyzing the more basic Hubbard-like Hamiltonian to 
show how it gives rise to the spin Hamiltonian. 

From the form of Eq. (JIJ, it is clear that the spin- 
phonon Hamiltonian we seek must be of the form 

V sp -p h = 2222b a lk S a (i)Sp(j)Q 7k , (2) 

where S(i) is the vector spin operator for site i and Q lk is 
the fcth normal mode amplitude at zero wavevector which 
transforms like the 7 component of a first rank tensor 
(vector) . We will discuss the normal modes in some detail 
in the next section. To implement the interaction of Eq. 
(J2J) it is convenient to classify both the normal modes and 
the spin components according to their transformation 
properties. This interaction represents a linear potential, 
i. e. a force on the phonon coordinate Q lk . Up to 
quadratic order in the displacements the terms in the 
elastic potential energy V D \ which depend on the Q lk are 

^ = 5E w X+Vph- (3) 

Minimization of this elastic energy with respect to the 
phonon coordinates leads to a phonon displacement 
which is proportional to the product of two spin func- 
tions, whose symmetry we analyze below. Furthermore, 
since the displaced ions carry an electric charge (albeit 
an effective charge), these displacements will give rise to 
a spontaneous polarization providing the necessary spin 
components are nonzero. This calculation will recover 
the symmetry properties of the phenomenological treat- 
ment. 



III. ZONE-CENTER PHONONS; NEUTRON 
SCATTERING MEASUREMENTS AND FIRST 
PRINCIPLES CALCULATIONS 

A. Generalized Displacements 

Since the normal modes are complicated linear combi- 
nation of atomic displacements we start by giving a qual- 
itative discussion of generalized displacements (GD's), 
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TABLE III: Irreducible representation of the paramagnetic 
space group of NVO. The vector representations are Bi u , B2u, 
and Biu which transform like z, y, and x, respectively. 



which are linear combinations of atomic displacements 
which transform according to the various irreducible rep- 
resentations (irreps), whose characters are given in Table 
IIIII for the paramagnetic space group of NVO. The ac- 
tual phonon modes are linear combinations (which we 
will give later) of these basis functions or CD's. 

We are only interested in CD's which transform ac- 
cording to the vector irreps B nu , for n = 1,2,3. We now 
give a qualitative discussion of these vector CD's. It is 
actually quite easy to generate the GD's which trans- 
form like vector components x, y, and z. First of all, one 
sees that assigning all atoms of a given Wykoff orbit the 
same displacement along the a-axis will give a GD which 
transforms under the symmetry operations of the space 
group (given in Table like the a component of a vec- 
tor. Since NVO has 6 crystallographically inequivalent 
sites this construction gives six GD's along each of the 
three coordinate axes, which we denote x n , y n , and z n , 
for n = 1, 2, ... 6. We now discuss the construction of the 
less trivial GD's which transform like a vector and which 
are shown in Figs. [21 EI and0] To generate these results, 
one can start with an arbitrarily chosen site to which 
a vector displacement along one of the three coordinate 
axes is specified. Then one generates the displacements 
of the other sites in the Wykoff orbit so as to reproduce 
the desired transformation properties. For example, to 
construct a z-like mode on the spines, we could assign 
the lower left spine site (in the lower left panel of Fig. 
|2J) a displacement along the x-axis. To be a z-like mode 
the pattern of displacements should be even under m Xl 
which fixes the displacement of the lower right spine to 
be that shown. Such a z-\ike mode should be odd un- 
der a two-fold rotation about an x-axis passing through 
the center of the cell. Applying this operation to the two 
lower spine sites fixes the orientation of the displacements 
of the upper spine sites. The other four symmetry op- 
erations give these same displacements. Had we started 
with a spine site with a displacement along the y-axis, we 
would have gotten a null displacement because this sym- 
metry with displacements along the y-axis is not allowed. 
Had we fixed the first site to have a displacement along 
the z-axis we would have found the trivial mode in which 
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FIG. 2: Generalized displacements yj, 27, z$, and X7 which 
transform like the components of a vector, for cross-tie sites 
(upper panels) and spine sites (lower panels). Atomic dis- 
placements (for the GD's indicated by the labels) in the x-z 
plane are represented by arrows, whereas those in the +y di- 
rection (— y direction) are indicated by filled (open) circles. It 
is easy to see that these modes couple to the uniform displace- 
ments. For instance, consider the mode z% shown here and in 
particular consider how the ionic displacements of the spine 
sites which are shown affect the cross-tie sites (not shown) at 
the corners and centers of the square. Imagine the ion-ion 
interactions to be repulsive. Then the nearest neighbors at 
negative 2 relative to each cross-tie get closer to the cross tie 
and the nearest neighbors at positive z relative to each cross 
tie get farther away from the cross tie. Thus, all cross-ties are 
squeezed towards positive z. This same reasoning also shows 
that even though this motion is confined to the x-direction, 
it will induce a dipole moment along the 2-direction. 
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FIG. 3: As Fig. [2] Here we illustrate schematically the vector 
GD's 29, 210, 2n, ys, 2/9, and j/10 for f sites. The placement 
of the sites reproduces the symmetry of an f site and is not 
quantitative for either V or O atoms. Three distinct f sites 
are occupied, one by a V atom and the other two by O atoms. 



all sites are displaced in parallel. The other GD's shown 
in the figures were generated in this way. Of course, the 
actual normal modes (phonons) consist of linear combi- 
nations of basis functions having the same irrep label and 
they will be discussed in the next subsection. 




In this section we present inelastic neutron scatter- 
ing (INS) measurements of the phonon density of states 
(DOS) in NVO along with the first-principles calcula- 
tions of the zone-center phonons. We will identify the 
phonons which have the correct symmetry to induce a 
spontaneous polarization. We will also attempt to esti- 
mate the local distortion which gives rise to the observed 
dipole moment in NVO. 

The INS measurements were performed using the fil- 
ter analyzer spectrometer (FANS) located on beamline 
BT4 at the NIST Center for Neutron Research^ . For 
energies above 40 meV, a Cu(220) monochromator, sur- 
rounded by 60' — 40' horizontal collimation and combined 
with a cooled polycrystalline beryllium filter analyzer 
was used. For the low energy spectrum (i.e. E < 40 
meV), a graphite (PG) monochromator with 20' — 20' 
collimation was used. The relative energy resolution of 
the FANS instrument is approximately 5% in the energy 
range probed. The powder Ni 3 V20s sample (about 20 
grams) was held at 12 K (paramagnetic phase) and 8 K 
(HTI phase) with a helium-filled aluminum can using a 
closed-cycle He 3 refrigerator. 

The first-principles total-energy and phonon calcula- 
tions were performed by the plane-wave implementation 
of the spin-polarized generalized gradient approximation 
(SP-GGA) to density functional theory (DFT)iSi We 
used 4x4x3 k-points according to Monkhorst-pack 
scheme and Vanderbilt ultra-soft pseudopotentials for 
which a cutoff energy of 400 eV was found to be enough 
for total energies to converge within 0.5 eV/atom. We 
considered the primitive unit cell of the NVO which con- 
tains 26 atoms as listed in Table II. Experimental lattice 
parameters were used in the calculations but the atomic 
positions were optimized to eliminate the forces down 
to 0.02 eV/A. The optimized positions are listed in Ta- 
ble. II, which are in excellent agreement with the experi- 
mental positions. Using the optimized structure, we next 



FIG. 5: The observed and calculated INS spectrum. The 
black vertical bars show the B2 U phonons whose intensity is 
proportional to the induced dipole moment when the system 
is distorted by the zero-point rms values of the modes. 



calculated the zone-center phonons and the correspond- 
ing INS one-phonon spectrum as described in Ref. \2'A 

The measured INS spectrum along with the calculation 
is shown in Fig.0 We observed almost identical spectra 
in the paramagnetic (T = 12 K) and the HTI (T = 8 K) 
phases and therefore show only the combined spectrum. 
The agreement of the calculations to the observed spec- 
trum is quite good, giving further confidence that the 
first-principle calculations capture the main physics. It 
also suggests that the phonon modes in NVO have small 
dispersion with wavevector. This is because the INS spec- 
trum is approximately averaged over a large range of Q 
and the calculations are only for Q = 0. The biggest 
difference between the INS spectrum and calculation is 
for the observed feature near 80 meV, which is calculated 
to be around 70 meV. As we discuss in detail below, in- 
terestingly this phonon has the right symmetry and the 
atomic displacement vector to induce a large dipole mo- 
ment. Hence, maybe the disagreement for the energy of 
this mode could be some indication of strong spin-phonon 
coupling. 

In order to identify the right phonon modes that can 
induce the observed dipole moment along the b-axis in 
NVO, we carried out the symmetry analysis of the zone 
center phonons. Table ITTT1 shows the character table for 
the irreducible representations (irreps) of the group G v 
for optical phonons at zero wavevector v. (For a review 
of group theory see Ref. 123 .) There are 26 atoms in the 
primitive unit cell and the representation T u induced by 
the vector space of these 26x3 = 78 atomic displacements 
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has the decomposition 

T u = 10A g + 8A U + 8Bi g + 13Bi u 

+7B 2g + 12B 2u + llB 3g + 9B 3u . (4) 

One can check that the vector representations which 
transform like x, y, and z are B 3u , B 2u , and B\ Ul re- 
spectively. To discuss the spontaneous polarization these 
arc the only irreps we will need to consider. Among 
the 78 phonons, twelve have ^> 2u symmetry, and can 
therefore produce the observed- spontaneous polariza- 
tion along the b-axis. One of the these twelve modes is 
acoustic (i.e. all atoms move uniformly along the b-axis) 
and will not be considered any further. To calculate the 
phonon energies and wavefunctions we found the eigen- 
values lj 2 of the matrix W(g = 0), which is related to 
the Fourier transform of the potential energy matrix by 
W T , arr ',af(q = 0) = VM^V T ,a; T ',a'(q = 0)y/M^, where T 
labels sites within the primitive unit cell, a labels Carte- 
sian components, and M T is the mass of the atom at site 
r. The corresponding eigenvectors oi™a are given in Ta- 
ble [W] Thus the displacements within the unit cell are 
given in terms of the normal modes operators Q n as 

Ur,a = J2°rt M T 1/2 Qn 
n 

= E°S\/^(4 + «n), (5) 

^ V 2M T Lu n 

where a) n is a phonon creation operator. The energies of 
the eleven optical modes are shown in Fig.0by black bars 
whose height for mode n is proportional to the average 
polarization Py™- 1 of that mode along the y axis. The 
polarization vector of the nth mode is estimated from 
the following formula: 

P (n) = — Vo (n) M -1 / 2 (6) 

T 

where f2 uc is the volume of the unit cell and Q r ms = 
yjTi/2uj is the average zero-point fluctuation. As we can 
see from Fig. \5\ half of the H 2u modes induce relatively 
small dipole moments. This is due to the fact that for 
these phonons, atoms mainly oscillate along the c-axis 
and the ^-component of the motion is only a second order 
effect. However for the other half, the motion is directly 
along the b-axis and therefore the induced dipole mo- 
ment is significant. Animations of these modes and more 
information can be obtained at Ref. We note that 
two particular phonons, one at 36 meV and the other 
around 70 meV, induce a significantly large dipole mo- 
ment. The eigenvectors of these and the other B 2u modes 
along with the magnitude of the induced dipole moments 
are given in Table llVl 

Figure [5] shows schematically how the oxygen atoms 
move in these two particularly interesting phonons. For 
the low-energy mode at 36 meV, the two oxygen atoms 
connecting the spine-spins move in the same direction. 




FIG. 6: A schematic representation of the top view of two 
particularly interesting B2 U modes whose displacement vec- 
tors are given in Table ITVI 



Therefore while one of the Ni-O-Ni bond angle decreases, 
the other Ni-O-Ni bond angle increases. Hence at first 
order, we do not expect large changes in the Ni-O-Ni 
superexchange due to this phonon. On the other hand, 
for the E = 70 meV mode, only one oxygen (which is 
connected to the cross-tie Ni spin) moves along the fa- 
axis. Hence, in this case, only one of the Ni-O-Ni bond 
angles changes from nearly 90° and therefore we expect 
this phonon to have important effects on the Ni-O-Ni 
superexchange interaction. Interestingly, the biggest dis- 
agreement between the experimental data and the cal- 
culated phonon energies happens for this phonon, which 
further suggests that it may have strong spin-phonon cou- 
pling. 

Finally, we discuss what kind of distortion is needed 
in order to induce the experimentally observed dipole 
moment whose magnitude is about P e xp = 1-25 x 10~ 4 
C/m 2 . We note that this induced dipole moment is 
much smaller than the calculated rms dipole moment 
(P rms - 46 - 70 x 10" 4 C/m 2 ) listed in Table Hvl This 
suggests that the local distortion should be at the or- 
der of P e xp /-Prms, which suggests that atoms move less 
than 0.001 A. This is a quite small distortion and would 
be very difficult to observe directly by neutron powder 
diffraction. 
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Ni(l) 
q = 0.90e 


Ni(2) 
q = 0.86e 


V(l) 
q = 1.18e 


0(1) 
q = -0.62e 


0(2) 
q = -0.68e 


0(3) 
q = -0.59e 


(meV) 


Wrms 

(A) 


(lO" 4 C/m 2 ) 


Mode 


(0,2/, z) 

y z 


(0,1/, 0) 

y 


(0,2/,z) 

y z 


(0,2/, z) 

y z 


(0,y,z) 
y z 


(x,y,z) 

x y z 








#4 


0.023 -0.446 


-0.119 


0.040 0.169 


-0.032 -0.015 


0.025 -0.314 


-0.024 0.005 0.048 


9.2 


0.47 


2.5 


#16 


0.551 -0.007 


-0.209 


-0.094 -0.063 


-0.099 -0.111 


-0.061 0.069 


0.013 0.073 -0.001 


21.3 


0.31 


1.7 


#27 


-0.205 0.058 


-0.110 


0.206 0.172 


-0.006 -0.354 


0.023 -0.061 


0.071 0.041 0.009 


30.4 


0.26 


3.3 


#29 


-0.008 -0.428 


0.019 


0.010 -0.265 


0.090 0.050 


-0.053 0.086 


-0.020 -0.040 -0.178 


31.1 


0.26 


6.1 


#34 


-0.164 0.066 


-0.286 


0.073 0.049 


0.109 0.163 


0.103 0.077 


-0.073 0.188 -0.068 


36.2 


0.24 


66.1 


#40 


0.178 0.113 


-0.166 


0.162 -0.012 


0.157 0.076 


0.171 -0.012 


0.029 -0.237 0.007 


38.9 


0.23 


18.7 


#49 


0.069 0.307 


0.021 


-0.056 -0.138 


0.023 0.057 


-0.069 -0.348 


-0.068 0.015 -0.139 


44.4 


0.22 


0.5 


#53 


-0.014 -0.022 


-0.037 


-0.043 0.165 


0.143 0.182 


-0.112 -0.080 


0.212 0.066 0.155 


49.8 


0.20 


16.3 


#64 


0.037 0.041 


-0.030 


0.264 0.110 


0.015 0.047 


-0.396 0.042 


-0.038 -0.027 -0.010 


69.2 


0.17 


46.1 


#70 


0.011 -0.011 


0.001 


-0.155 0.017 


0.365 -0.169 


-0.081 0.011 


-0.143 -0.004 0.090 


93.3 


0.15 


23.0 


#78 


0.010 0.019 


0.004 


-0.089 0.238 


0.150 -0.099 


-0.029 0.000 


0.203 0.009 -0.187 


103.0 


0.14 


11.3 



TABLE IV: Mass-weighted atomic displacements of B2u phonons which induce a dipole moment along the b axis (normalized 
so that the sum of the squares of the components equals unity for each mode). (The acoustic B2 U is not tabulated.) Each 
component of the mass-weighted displacement represents the atomic displacement divided by the square root of the respective 
atomic mass. The mass- weighted displacements are given for the sites listed in column #2 of Table ITT1 The displacements of 
the remaining atoms in each Wyckoff orbit are fixed so that the mode transforms like B2 U , i. e. like the y-component of a 
vector. The calculated atomic charges, magnitude of the rms displacement Q r ms (where m p is the proton mass), the rms dipole 
moment, and the mode energy is also given. 





i 


2,x 


m y 


rh z 


r 1 


i 


i 


1 


1 


r 2 


l 


i 


-1 


-1 


r 3 


i 


-i 


1 


-1 


r 4 


i 


-i 


-1 


1 



TABLE V: Irreducible representations of the group G v for the 
incommensurate magnetic structure with v = (q, 0, 0). Here 
it is simplest to use the symmetry operations fh y and rh z , 
such that rh y v — (x,y + ^, z + ^) and m z r — (x, j/+ §,2+ |)- 

IV. SYMMETRY OF THE SPIN STRUCTURE 

In Refs. 0, 0, 0, an d HU the application of represen- 
tation theory to the determination and characterization 
of magnetic structures is discussed in detail. In Table 
we give the character table for the irreps for an arbitrary 
wavevector of the form (q, 0,0). Because the irreps are 
one dimensional, each spin basis function is an eigenvec- 
tor of the symmetry operator with the listed eigenvalue. 

In these references it is shown that the HTI phase is 
described by a set of five complex amplitudes associated 
with the irrep T4. Here we call these a s x , ia s>y , and a SjZ 
to describe the orientation of the spine spins and a c<y 
and a Cj 2 to describe the orientation of the cross-tie spins. 
When the LTI phase is entered additional variables as- 
sociated with the irrep T2 become nonzero. These LTI 
variables are here denoted ib s ^ x , b s ^ y , and ib StZ to describe 
the orientation of the spine spins and 6 C Z to describe the 
orientation of the cross tie spins. Because the crystal is 



centrosymmetric, it is showrA£i2£ that within a given rep- 
resentation all these complex structural parameters (the 
a's and &'s) can be written in terms of a real amplitude 
times a complex phase factor which is the same for all 
variables of the same irrep, I^ or T2, in the sense that 

at, a = ±\a t . a \e l ^ TI , bt, a =±\b t . a \e^ LTI , (7) 

where t denotes either spine or cross-tie. It is further 
expected 4 *^ that </>hti — </>lti is ±7r/2. Thus in these 
two phases we may use the results of Table VIII in Ref. 
4 to write the spin components of the six Ni ions in the 
unit cell as 

fiW(Ri) = (<».., • ^. , ),'<■ " 

+«* - <'<;,)< ,<, ,( (8) 
+(-Kv + & :, a K iq,Ri ( 9 ) 

S«(Ri) = • ,7,... " 

(10) 

si 2 Hn 2 ) = '/..,. • '7,.., • 

+(-<x-«C>- <q - Ra (11) 

4 2 >(R 2 ) = " ' 

+(-Kv ~ ^K"^ 2 (12) 
5p)(R 2 ) = (a s , z - ib s , z )e^ 

+ « z + z&;je-^ (13) 
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(14) 


4 3) (R 3 ) 








+«, + ^)^ iq ' R3 


(15) 


5i 3 )(R 3 ) 










(16) 


&W(R 4 ) 


( ^s.x ~~ ^s,x)^ ^ 






+(-<x + ^,Je- iq ' R4 


(17) 


4 4) (R 4 ) 










(18) 


S( 4 >(R 4 ) 








+« 2 -^:>- iq - R4 


(19) 


Si 5 \R 5 ) 




(20) 


4 5) (R 5 ) 


= a^e^+a^e"^ 


(21) 


5l 5 )(R 5 ) 


= + < z e- iq R5 


(22) 


5i 6 )(R 6 ) 


= -6 Cia; e lq - R6 - 6^e- iq - R6 


(23) 




iq-R.6 ~* ~— iq*R-6 


(24) 


^ 6) (R 6 ) 


= a c , 2 e q + a e 4 


(25) 



Here the superscripts on 1, 2, 3, 4 on S label the spine 
sites and 5 and 6 label the cross-tie sites c and c', respec- 
tively, (as in Fig. [7] below) and R n is the position of a 
spin n. Also q = qi with q ss 0.28(27r/a). 2,4 



V. RELATIONS FOR THE STRAIN 
DERIVATIVE OF THE EXCHANGE TENSOR 

In this section we obtain explicit forms for the most 
important spin-phonon coupling matrices. For this pur- 
pose we start by introducing notation for the principal 
exchange interactions. We write the interactions between 
spins on sites i and j as 



= ^2x al3 (i,j)S a (i)Sp(j) , 

a/3 



(26) 



where X a p(i,j) — Xp a {i, i), of course. For nearest 
neighbor (nn) interactions between spine spins we set 
Xafiihj) — U a p(i,j), for next nearest neighbor (nnn) 
interactions between spine spins we set X a p(i,i) — 
V a fj («, j) , and for nn interactions between spine and cross- 
tie spins we set X a/ s(i,j) = W a fj(i,j). We may further 
decompose the exchange tensor into its symmetric and 
antisymmetric parts. For example, for nn spine-spine in- 
teractions we write (omitting the site labels i and j) 



U 



J xx J xy H~ D z Jxz 



J, 

■3 X 



J y 



■Jy 



Dy 



Dy Jy Z 



Jz 



(27) 



where D is the Dzyaloshinskii-Moriya (DM) vector&S 
Similar decompositions will be made for V and W in 



terms of symmetric tensors K and L, respectively, and 
the DM vectors E and F, respectively. 

Now we consider the gradient expansion of these ex- 
change tensors. From the above development we write 
the displacement u of the rth atom in the unit cell in 
terms of the amplitude Q lp of the GD labeled -f p as 



i Q (R,r;7 p ) = Q lp a(~f p ;T,a) , 



(28) 



where a(7 P ; r, a) is the a-component of the displacement 
of atom r in unit cell at R for the GD labeled j p (which 
are shown in Figs. lUGl andQJ. If Z represents a com- 
ponent of an exchange tensor, then we write 



dZ 



1p 



E 



dZ 



a{l p -,T, a) 



(29) 



The aim of the present paper is to determine which such 
derivatives are required to completely determine the tri- 
linear spin-phonon coupling. The actual calculation of 
these coefficients will be undertaken separately. Thus, 
for each GD 7 P , we consider the interaction 



— Q 



-fe 2^ 2^ — — s a (i)s {j) . 

a (3 ij 



dQ 



(30) 



7p 



Our objective is to express the results for the spontaneous 
polarization due to the trilinear coupling in terms of the 
parameters dX a p{i 1 j)/dQ^ p . 

Here we analyze the gradients of the nn interactions 
between spine spins. (Similar analyses of next-nearest 
neighbor spine interactions and of nn spine-cross tie in- 
teractions are given in Appendices.) We introduce the 
coupling between spine sites #1 and #4 in Fig. |7| 



dU aP {lA) 



dQ 



1p 



u af3 



(31) 



Clearly U 7p has to be invariant under the symmetry op- 
erations of the crystal. But m x takes the bond in question 
into itself (and interchanges indices). This indicates that 
these interaction matrices must satisfy 



(T x \J Xp cr x = -\J Xp 
(T x U yp (T x = t>» 

(T X U Z "(T X = LP* , 

where tilde indicates transpose, 



-10 
1 
1 



(32) 



and later 



1 
0-10 
1 
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y = 
y = d 



3' 






— J x 

u xy 


— J x 

J xz 


— J x 

J xy 





D x 

x 


— J x 
°xz 


_ D x 

x 






(38) 



y = 



y = -d 



y = 



4' 



O 

1 



FIG. 7: Diagram of an a — c plane used to specify nearest- 
neighbor and next-nearest-neighbor interactions along a single 
spine. Circles are spine sites and square are cross-tie sites and 
d = 0.136 (see Table HTl . The dashed rectangle indicates the 
unit cell. Interactions in other a — c planes are obtained by 
using translation symmetry. 



dU a p(2,3) 



dQ 



Up 



dQz p 



-a-yCr z JJ Vp cry<T z 



iJ .y. rf D^. 



- T v 

yy 

-Jjjz 



XX 

-jyy 

Dv 



-D v 

v 

- l y 

~J V ~ 



-J. 



XX 



D: 



D z 

^ z 

- l z 

Jy 



(Ty(T Z lJ Z "<Ty<T Z 



-D z 
v 

- J z 

• J yz 
~ J zz 



(39) 



(40) 



1 
1 

0-1 



(33) 



Here we used the fact that the GD's were constructed to 
have known symmetry: 



and 



m a Qt3 p = (1 - 25a,p)Qf3 p 
2 a Qi3 p = (-1 + 25 a .p)Qp p 



(34) 
(35) 



where 8 a ,b is unity if a = b and is zero otherwise. 
In view of Eq. I1-S2H , we have 








JX 

u xy 


JX 

,J xz 


LP" = 


JX 

J xy 





D x 

X 




jx 

^ xz 


~Dx 







jy 

^ XX 


DV 


—D y 

y 


W" = 


-Dl 


TV 

yy 


jy 

J yz 






Jyz 


Jzz 




J Z 

w XX 


D z 

z 


-D z 
y 


U z » = 




J z 

yy 


T z 

J yz 




. D y 


T z 

J yz 


l z 

zz 



(36) 



where the index p on the superscripts of J and D are 
left implicit and and DZ P (and similarly later for 
superscripts on K, L, E, and F) are defined to be 

dJ a p(lA) 



t7j> 
J af3 



DZ^ dD A 1A) - (37) 



dQ lp ' a ' dQ lp 

Then we obtain the #2-#3 interaction from the above 
by 2 X , a two- fold rotation about the x-axis, so that 



dU afi (2,2) 
dQx v 



= (Ty(T Z U Xp (Ty(T Z 



where we used Eq. (|35p. 

We obtain the #4-#l' interactions by applying the 
glide operation m y to the #2-#3 interaction, so that 



dQx v 



dU(2,3) 
' dQ Xv ' 






JX 

u xy 


— J x 

• J xz 


JX 





_ D x 

^ X 


— J x 

°xz 


D x 

^ X 






(41) 



dU a p(4, 1') 



dQ 



Up 



<9U(2,3) 
' J y xx Dv Dv 

-d\ jy y -jy yz 



-D v 

y 



'Jyz 



(42) 



dU a p(4,l') 
dQz v 



<9U(2,3) 
cr — tt^ cr,, 



dQz 
- T z 

" XX 

D z 
^ z 

Df, 



-D z -D z 

^ z ^ y 

-JL J, 



J:, 



yz 
T z 



(43) 



and finally we get the #3-#2' interaction by applying a 
two-fold rotation about the x-axis to the #4-#l' inter- 
action to get 



dU a p(3,2') 
dQx p 



9U(4, 1') 
^"^QxT^" 



— T x 

u xy 
jx 

,J XZ 





D x 



XZ 

D x 



(44) 
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dU af3 (3,2') 



dQ 



dU(4, 1') 

-<T,,(T 2 — <X„(X , 



— j y 



9Qy P 



U z J yy 



-D y 

V 



Jyz 



(45) 



and 



dQ Zv 



au(4, i') 

-<J „<T , — — er „er , 



XX 

D z 

z 



dQ z 
-Dl 



Df 



T z —J 

Jyy J 



(46) 



VI. MEAN FIELD SPIN-PHONON 
HAMILTONIAN 

A. Mean Field Results 

Here we treat the nn spine-spine interactions in detail. 
Analogous calculations for the nnn spine-spine and nn 
spine-cross tie interactions are treated in an Appendix. 
We evaluate the spin-phonon Hamiltonian 7i 7p of Eq. 
(|30fl at the mean-field level. In other words, for the spin 
operators we simply substitute their average values as 
given in Eqs. i|8Jl-l|25fl. One sees that Tt Xp for nn spine- 
spine interactions, for instance, will consist of contribu- 
tions proportional to J X y, to J X S, and to D x p . To illus- 
trate the calculation we will explicitly evaluate the first 
of these, which we denote Ti. Xp (Jx y )' 

Hx p (Jxy) — Qx p J X y ^ ] ' ^ ^ ^x ' 



_ c(2) c(3) _ o(3) c(2) , c(4) a(l') 
^ x y x y i x y 

+ e(l')o(4) _ c(3)e(2') - S*( 2 ')s*( 3 ) 



(,47) 



where, since we included all interactions within a single 
unit cell, the sum is over all -/V uc unit cells. (In this 
summation only terms involving both q and — q survive.) 
Thus 



T~(-x p (Jxy) — 2N uc Q Xp J x p e 



-iqa/2 



x [(a SlX + ib StX )(ia* s 



"s,yj 

+ b s , y )[-a* +ib*) 



f ib s<x )[ia* y + b* y ] 

b s,y){alx +iK,x)] 



-(-as, 

-{ia s , y 
= 16N uc Q Xp J* p cos{qa/2) 
x %[a* StX a StV + b s<x b*J . 

The other terms proportional to Q Xp are 

n Xp (J xz ) = lW uc Q Xp J*i sm(qa/2) 
x3[a StX a* + b StX b* ] . 



(48) 



(49) 



H Xp {D x ) = WN uc Q Xp D%» cos{qa/2) 
x%{-a Sty a* +b s>y b* ] 



(50) 



In view of Eq. all these terms involving Q x vanish, as 
was found from the phenomcnological formulation. Sim- 
ilarly, all the terms involving Q z also vanish, again in 
conformity with the phenomenological argument. 

We are thus only left with terms involving Q y . In 
Appendix C we find that the strain dependence of the 
nn spine interactions give 

Hy, = 16N uc Q yp J2 A i" n) 3 K A"] , (51) 



fi,f— X ,y,z 



where 



J Vp r 


D v z p s 


Dy P C ' 


-D y /s 


- 7 Vp r 

JyyC 


- l Vp H 
,J yz o 


D V y"C 


Jy z S 


-Jllc 



(52) 



where c = cos(qa/2) and s = sin(qa/2). Using the results 
of Appendix C we find that the nnn interactions give a 
result of the form of Eq. (|51|l but with 



K xx c 
E v z p s' 
-K y x p c' 



-E, 



Vp s' 



E 



Vp r l 

yy c 
Vp j 



-El p s' 
-Klld 



, (53) 



where c' = cos(ga) and s' = sin(qa). Using the results 
of Appendix C we have the results for the spine-cross tie 
exchange gradients: 

V Vp = 16N uc Q yp I E A a^M^] 

\a=y,z f3 



ct—x (3 



(54) 



where A sx is 



L xx s 



T u - a" 

J-lyyb 



[Ll p v - F v z p ] c" 
[Lll + F^] c" [L v y l-F^] 



[L% + F y z p ] c" [Lll - F y y p ] s" 
[L v v p z + F yp ] c" 



L v z p z c" 



(55) 



where c" = cos(ga/4) and s" = sin(ga/4). These results 
agree with the phenomenological model, in that V Vp is 
only nonzero when both the "a" and the "6" irreps are 
simultaneously present and they can not have the same 
phase (lest a* b be real). 



B. Summary 

Here we show how the above results lead to an evalu- 
ation of the spontaneous polarization. If we combine the 
results of Eqs. l(5TJl . (|53^l . and (|54*jl . we may write 



V, 



Up 



NucQy p ^y p i 



V 



1,12 



(56) 



11 



But the GD's are related to the normal modes via 

Qv p = E o^M-^Q n , (57) 

n=l,ll 



+Z 



-z 



where the orthogonal transformation o is essentially that 
given in Table IIVI apart from normalization factors due 
to replacing a Wyckoff orbit by one of its atoms. Note 
that the GD y p is associated with atoms all having the 
same mass M p . The magnetoelectric trilinear interaction 
in terms of normal mode operators is thus 

V = N^Ki^Mp^Qn , (58) 

n,p 

Since the unperturbed elastic energy is Ha — 
(N uc /2)J2 n UJ nQn: we see that the trilinear interaction 
induces the displacements 



Wn)=EVtt" 1/2w n 2 , 



(59) 



from which we obtain the induced spontaneous polariza- 
tion as 



Py = 



« 'lie 



(60) 



where the atomic charges q p are given in Table llVl The 
only ingredients we do not have for this evaluation are 
the various gradients of the exchange tensors. 



VII. TOY MODELS 

In this section we will consider two toy models. In the 
first one, we consider a single spine with one Ni per unit 
cell, but with two oxygen atoms symmetrically placed 
on either side of the spine. In this version, all atoms lie 
in a single plane. In the second version, the size of the 
unit cell is doubled. In one plaquette the oxygen atoms 
are both displaced equally perpendicularly to the original 
atomic plane and in the next plaquette the oxygens are 
oppositely displaced. In this version one of the mirror 
planes becomes a glide plane. These models illustrate the 
simplifications which arise when the system has higher 
symmetry than the buckled kagome lattice of NVO . 




FIG. 8: The two toy models. In each case the unit cell is 
bounded by a dashed line. The filled circles are Ni sites and 
the squares are O sites. 



where = dJ a p/dQ lp , D1 = dD a /dQ lv , and the in- 
dex p is left implicit. The Hamiltonian is invariant under 
mirror reflections with respect to each coordinate axis. 
Taking account of the symmetry of the displacement co- 
ordinate and the fact that m x interchanges indices of the 
exchange tensor, we have that 



<T a .J 7 (T 
(Tj ) J 7 (T 
(7 - .7 (7 



(1 - 25 X)7 )JT 
(1 - 28 Vtl )r 
(1 - 2£ Z , 7 )JT 



(62) 



As a result of this symmetry we have that JP 

j» = 



D\ 
-Df 




—D 

D z y 



(63) 



We construct the trilinear spin-phonon interaction, V, as 

V = J2(Qy P DV z PC y + Qz P D y PC *) > ( 64 ) 



where 



C a = Y^[S x (n)S a {n + l)-S a (n)S x {n + l)]. (65) 



A. Unit cell with two atoms 

In this section we consider the toy model shown in the 
left panel of Fig. |H1 We first analyze the symmetry of 
the strain derivatives of the exchange tensor. By trans- 
lational symmetry all nearest neighbor interactions are 
equivalent. So we set 

dJ a p{n, n + 1) 



dQ 



1p 



_ t7p 
— j q/3 



XX 



Jxy + ^xz 



Jl 



yy 



Jyz 



D" 1 

V 

DJ 



Jxz + Dy Jyz 



J]z 



(61) 



Now we replace the spins by their equilibrium values. 
We note that each spin component belongs to a separate 
representation and we write 



S x (n) - S x (qy nqa + S x (q)*e- in " a , 
S v (n) = S y (q)e in * a + S v (q)*e- in * a . 
S z (n) = S z {qy ma + S z (q)*e- in « a 



(66) 



Then, keeping only those terms which survive the sum 
over n we have 

V = 4Nsm{qa) ^ [Qy p D v /r x (q)r y (q) sin(<p x - <f> y ) 
p 

+Q Zp D*»r z {q)r x {q) sin(0 x - ^)] , (67) 
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TABLE VI: Basis spin functions for sites #1 and #2 in the 
unit cell in terms of the complex-valued Fourier components 
S a (q) for irreps characterized by the eigenvalues of m y and 
the glide operation m z . 



r 


m v 


m z 


s(#i) 


S(#2) 




+ 


+ 


(S x (q),0,S z (q)) 


(&(<?), 0, "&(</)) 


r 2 


+ 




(S x (q),0,S z (q)) 


(-S x (q),0,~S z (q)) 


r 3 




+ 


(0,S y (q),0) 


(0,S y (q),0) 


r 4 






(0,S y (q),0) 


(Q,-S v (q),0) 



where N is the total number of Ni spins and we set 
S a (q) — r a (q)e l ' l>0 ' , where r a (q) is real. As found 
befor o 3 i i this interaction is only nonzero when a) two 
different representations are condensed and b) the two 
representation have different phases 4>. In view of our 
previous results, it is clear that the appearance of only 
strain derivatives of the DM vector is an artifact of the 
rather high symmetry of this coplanar model. 



B. Unit cell with four atoms 

Now we consider the noncoplanar toy model shown 
in the right panel of Fig. |S| The first two symmetry 
relations of Eq. (|62[1 remain valid, but the third one now 
results from the glide plane which involves a translation 
along the chain. If J" denotes the strain derivative of 
the exchange tensor for coupling sites 2n and In + 1, and 
J"? that for sites 2n + 1 and 2n + 2, then we have 



±J* 




±J% Z 



(68) 















J a ± = 


-Dl 







(69) 





















-D z 

V 




Ji = 





+ T* 

yy 





(70) 




D z 

. V 










Next we characterize the spin structure. There are four 
irreps for which the basis functions are listed in Table IVTl 
Thus we write for spins #1 (Si) and #2 (S2) in the unit 
cell, 

S hx (X) = [SP(q) + Si 2 \q)}e^ x 

+ [^ 1 )(g)*+5( 2 )(g)*] e ^ x , (71) 
S 1>y (X) = [S^(q) + S^(q)}e i " x 

+ { S W( q y+S^( q y}e-^ x , (72) 




+z -z 



+z -z 




+z -z 




x +z 



+z 




FIG. 9: Generalized displacements which transform like vec- 
tors. Open circles (circles with inscribed "x") are displace- 
ments out of (into) the page. The sites at positive or negative 
z are indicated. The upper panels show modes in which the 
atoms move only in the z-direction. Upper left: an z-like 
mode X3. Upper right: a y-like mode j/3. The lower panels 
show z-like modes. Left: 23 with motion only along the x-axis 
and right: 24 with motion only along the j/-axis. 

s lt ,(x) = [sP(q) + si 2 H q )}e^ x 

+ [S«( 9 r +S( 2 >( g )*]e-"? X , (73) 
S 2 , X {X) = [S^{q)-S^{q)]e^ x 

+[Si 1 \qy-SM(qy]e- i « x , (74) 
S 2<y (X) = [S^(q)-S^(qW qX 

+ [S y 3) (qy-S^( q y]e-^ x , (75) 
S 2 AX) = [-S^{q) + Sf\q)}e^ x 

+[Si 1 \ q y +S^{qy\e^ x , (76) 

where the superscript labels the irrep as in Table l*VTft . 
Thereby we find the trilinear spin-phonon coupling (when 
the spin operators are replaced by their values: 



V = 4iVsin(ga)^ 



UQ Xp J2 + Q Vp (yj*> + WD*>) 



(77) 



where 



U = <z(s£Xqys£Xq) + S£Xq)S?\ q y) (78) 



V = ^[s^(q)Sf\qy + SP{q)S y 5 \qy) (79) 

[s^\ q )sf\ q y + sf\ q ys^{ q )) (so) 



W = 9 



13 



X = 3 



(sW( q ysl 2 H q ) + sP( q )sP( q y) (si) 

Y x = %(sP( q ysW( q )) (82) 
Y y = 3(4 3 )( g )*4 4 )( g )) 
9f(5W(«)SW(«)* 



(83) 



(84) 



The general symmetry arguments indicate that there can 
not be a polarization along i. We see that J7 van- 
ishes because all the components within a single rep- 
resentation have the same phase, so that, for instance, 
sP(q)si 1 \q)* is real. Here we see that, depending on 
the spin structure the spontaneous polarization can ei- 
ther be along j (if either both irreps #1 and #3 are 
active or both irreps #2 and #4 are active) or along k 
(if either both irreps #1 and #2 are active or both irreps 
#3 and #4 are active.) These results are exactly what 
the phenomenological analysis would give. 



VIII. CONCLUSIONS 

In this paper we present neutron scattering measure- 
ments of phonons in NVO, the first-principles computa- 
tion of the zone-center phonons and their symmetry anal- 
ysis. We identified two particularly interesting phonons 
among the twelve B2« modes which have the right sym- 
metry to induce the experimentally observed dipole mo- 
ment along the b-axis in NVO. Using the calculated 
atomic charges and the eigenvectors we conclude that the 
required distortion to induce the observed dipole moment 
is small ( 0.001 A) and would be difficult to observe di- 
rectly by neutron powder diffraction. Finally, we present 
a symmetry analysis to characterize the microscopic mag- 
netoelectric coupling in NI3V2O8. The method can easily 
be applied to similar systems such as TbMnC>3. The re- 
sult of this analysis is that we can specify those strain 
derivatives of the exchange tensor which should now be 
the targets of more fundamental quantum calculations, 
perhaps based on the LDAii or similar schemes. (These 
results are given in Eqs. H52|) - (|55ll . where the superscript 
indicates the component of the gradient.) In Subsec. 
IVI Bl we show how these gradients lead to an evaluation of 
the spontaneous polarization from first principles. In Sec. 
IVI 1 1 we also studied some structurally simpler toy mod- 
els. A general conclusion is that the local site symmetry 
in NVO is low enough that almost all strain derivatives 
of the exchange tensor are involved. It is true that the 
Dzyaloshinskii-Moriya interaction must be present (and, 
of course, it is allowed in such low symmetry systems). 
In the presence of such interactions ferroelectricity results 
from the strain dependence of all the components of the 
exchange tensor, even that of the isotropic exchange in- 
teraction. 
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APPENDIX A: NEXT-NEAREST NEIGHBOR 
SPINE INTERACTIONS 

In this section we consider next-nearest neighbor (nnn) 
interactions between spins on a given spine. Since only 
the gradients with respect to Q Vp are needed, we only 
consider those here. We set 



dJapQ.,1') 



- v Vp 



dQ 



a/3 



(Al) 



The operation 2 y takes this bond into itself with reversed 
indices. So, taking account of the transformation prop- 
erties of Q Vp , we require that 



2yV V * 



yyp 



(A2) 



where 2 y is the two-fold rotation operator. In terms of 
matrices, this relation is 



where 2. 



2 y V y »2 y = V y " 
Thus 



K y 

XX 


E y 


K y 

xz 


-E y 

J -'z 


K y 

yy 


E y 

.r 


K y 


-E y 

*^x 


Ji zz 



We have that 



dJ a0 {2,2') 



dQ 



-2 x V y *2 x 



Up 



so that 



dJ a0 {2,2') 



dQ Vp 
We have that 



-K y 

XX 

-E y 

K xz 



E y 

-K y 

yy 

E y 



K y 

XZ 

-E y 

X 

-K y 



dJ a p(3',3) 



dQ 



= -<ryV lp ay , 



Up 



so that 



dJ aP (3',3) 
dQ Vp 

We have that 



-K y 

XX 

-Ev 

~K y z 



E y 

-K y 
yy 

-E y 



-K y 
E y 

X 

-K y 



SJ Q/3 (4',4) 



dQ 



= cr T V>< 



(A3) 



(A4) 



(A5) 



(A6) 



(A7) 



(A8) 



(A9) 
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FIG. 10: As Fig. |7| Here we show the eight different nearest- 
neighbor spine-cross tie interactions, which we label A ... H. 
The arrow points from the first site index to the second site 
index. We also give the symmetry operation one has to apply 
to interaction A to get each of the other interactions. 



so that 



3J Q/3 (4',4) 



dQ 



Vp 



K xx El 
-K y 

XZ 



E y 



-K y 

^xz 

-E y 

*^x 

K y 

Jv zz 



(A10) 



APPENDIX B: SPINE CROSS-TIE 
INTERACTIONS 

In this section we analyze the spine-cross tie interac- 
tions, in Fig. 1101 again keeping only derivative with re- 
spect to Q Vp . We set the interaction of type A to be 



dJ a p{c, 1) 



dQ 



Up 



so that (when the index p is suppressed) 



W y = 
We have that 

so that 
<9J(c,4) 



L y 

^xx 

L y - F y 

xy z 

py i py py 

xz 1 y yz 



py + py py - F y 

xy 1 z xz y 



L y 

yy 



py 



L y + F y 

yz ' x 

L y 

^ zz 



(Bl) 



(B2) 



<9J(c,4) 
dQ Vp 



= <X T W>< 



(B3) 



dQ 



Up 



We have that 



L y 

XX 

.pv + py py 

xy 1 z yy 

_ tv _ py tv — py 

xz y yz x 



.py _ py -pv + py 

xy z xz 1 y 

py i py 

yz ' 1 x 

L y 



dJ(c,2) 
dQ Vp 



= -2 x W y *2 x , 



(B4) 
(B5) 



so that 

aj(c,2) = 

dQ Vp 
We have that 

so that 
<9J(c,3) 



-L y 

xa 

py _ py _py 

xy z yy 

py _|_ py -py + py 

xz 1 y yz 1 x 



py _|_ py py _ f « 
.py _ 

yz A a; 

zz 



dJ(c, 3) 

9Qy P 



= -1W"1 



(B6) 



(B7) 



We have that 

so that 
<9J(c',l) 



sex 

_£y _|_ py -py 

xy 1 2 

.£« _ jry -py _i_ i?y 

X2 y yz 1 a; 



_ pv -py _|_ fy 

- F y 

yz x 

-L y 



9J (^) - 2 W2 

dQy p " ' 



(B8) 



(B9) 



9Qy P 
We have that 

so that 
<9J(c',4) 



L y 



. tv — py T,y _ py 

™' z xz y 



.py _i_ py py 

xy 1 z yy 



^xz ' y yz ' x ^ zz 



dJ(c',4) 
dQy P 



= er 2 W>< 



(BIO) 



(Bll) 



9Qy P 
We have that 

so that 
d3(c',2) _ 



L y . 



py + py ^py + py 

xy 1 z xz 1 y 



py _ f y £y 

* yy 



.£3/ _fy 

y2 a: 



_£y _ + py py 

xz y yz 1 a: 22: 



c>J(c',2) 

dQyp 



= -2 Z W>2 Z , 



(B12) 



(B13) 



dQy p 

We have that 

so that 
<9J(c',3) 



-L y 



-L y ,-F y Li. -FX 



.py _i_ py _py py + py 

xy 1 2 yy ys 1 x 



py i py py _ f y 

a)2 1 y y^ 33 



9J(c',3) 



-<T,W><X, , 



(B14) 



(B15) 



dQ 



Up 



-py pv + f y _^y + py 

xx xy 1 2 xz 1 y 

^3/ _ py _py py + py 

xy z yy yz 1 x 

-py - py py _ jro 

xz y yz x 22 
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APPENDIX C: SPIN-PHONON INTERACTION 



Now we collect terms involving J^ z : 



1. nn spine interactions 

Now we collect terms involving J% : 



(1)S X (4) - S X (2)S X (3) 



+S x (4)S x (l') - S X (3)S X (2') 



— Q y J xx e 



(flrr +ib x )(-a* x +ib* x ) 



-(-a x +ib x )(a* x +ib* x ) 
+{-a x - ib x )(a* x - ib* x ) 

~(«x - ib x ){-a* x - ib* x ) 



+ c. c. 



= 4Q V J V XX N UC (e^"K&* - iM*] + c. c.) 
= -16iV uc Q w J* x cos(ga/2)9f[o x 6*] . (CI) 

Now we collect terms involving M : 



Vy(Jyy) — QyJyy*^2 



Sy(l)Sy(4) - Sy(2)Sy(3) 



+ Sy(4)Sy(l')-Sy(3)Sy(2') 



QyJ V yyE< 



[ia y + b y )(ia* y - b* y ) 



-(ia y -b y )(ia* y + b* y ) 
+{-ia y - b y )(-ia* y + b* y ) 

-(-ia v + b y )(-ia* -b*) 



+ c. c. 



- AQ y jy y N uc (e iqx ^[-ia y b* y + ib y a* y ] + c. c.) 
= l6N uc Qyjy y cos(qa/2)<Z[ayb;} . ' (C2) 

Now we collect terms involving J| z : 



Vy(JL) = QyJLj2 



S z (1)S Z (4) -S z (2)5,(3) 



+S Z (4)S Z (V) - S Z (3)S Z (2') 



QyJ V zz J2< 



(a z + ib z ){a* z - ib* z ) 



-(a z -ib z )(a* z +ib* z ) 

+(a z + tb z )( a ;-tb;) 

-(a z - ib z )(a* z + ib* z ) 



c. c. 



- 4Q v J y zz N uc (e igXli [-ia z b* z + ib z a* z ] + c. c.) 
= 16N uc Q y jy z cos(qa/2)Q[a z b* z }. (C3) 



Vy(Jy Z ) — QyJy Z 



S v (l)S z (4) + S z (l)S y (4) 



-S z (2)S y (3) - S y (2)S z (3) 
-S y (4)S z (l>)-S z (4)S y (l>) 

+S y (3)S z (2') + S Z (3)Sy(2') 



QyJ V y Z J2' 



(ia y + b y )(a z - ib* z ) 



+ (a z + ib z )(ia* y - b* y ) 
-(a z - ib z ){ia* y + b* y ) 
-{ia y -b y ){a* z +ib* z ) 
-(-ia y -b y )(a* z -ib* z ) 
-{a z + ib z ){-ia* y + b* y ) 
+(—ia y + b y )(a* z + ib* z ) 



+ c. c. 



+ (a z - ib z )(-ia* y - b* y ) 

= AQ y jy z N uc (e^[a y b* z + b y a* z - a z b* y - b z a* y 
+c. c.) 

= -16N uc Q y S m(qa/2)jy z Q[a* y b z + b* y a 
Now we collect terms proportional to D v z : 



(C4) 



v y {Di) = Q y DyJ2 



S X (l)Sy(4) -Sy(l)S X (4) 



+S x (2)S y (3) - S y (2)S x (3) 

+ S X (4)Sy(l') - Sy(4)S X (l') 

+S x (3)S y (2') - Sy(3)S x (2') 



= QyD y z J2 e 



zqxi4 



(a x +ib x )(ia y - 6*) 



-(ia y + b y ){-a x + ib%) 
+(-a x + ib x )(ia* y + b* y ) 
-(ia y -b y )(a* x + ib* x ) 
+{-a x - ib x )(-ia* y + b*) 
-(-ia y - b y )(a* x - ib* x ) 
+ {a x - ib x )(-ia* y - b* y ) 



+ c. c. 



-(-ia v + b v ){-a* x -ib* x ) 

= AQyD v z N uc (e^ [-a x b* y - b x a* y + a y b* + b y a* x ] 
+c. c.) 

= -WN UC Q V sm{qa/2)D y z Q[a*yb x + b* y a. 
Now we collect terms proportional to D y : 



(C5) 



V y {Dl) = QyDlY, 



-S X (1)S Z (4) + S Z (1)S X (4) 



-S X (2)S Z (3) + S Z (2)S X (3) 
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+5,(4)5,(1') -5,(4)5,(1') 
+5,(3)5,(2') - 5,(3)5,(2') 



= QyD v y J2 e 



iqxi4 



-{a x +ib x )(a* z - ib* z ) 



+(a z + ib z )(-a* x +ib* x ) 
~(-a x + ib x )(a* z + ib* z ) 
+ (a z - ib z )(a* x + ib* x ) 
+(-a x - ib x )(a* z - ib* z ) 
-{a z + ib z )(a* - ib* x ) 
+{a x -ib x )(a* z + ib* z ) 

-{a z -ib z ){-a* x -ib* x 
+c. c. 

= 4Q„£){(JV UC {e^H\a x b z - b x a* z + a z b* x - b z a* x ] 
+c. c.) 

= -WN uc Q v cos(qa/2)Dp[a x b* z + b x a z } . (C6) 

2. nnn spine interactions 

First we consider terms proportional to K% x (where we 
omit the p index): 

Vy{Kl x ) = Q y KV x Y,[S x {X)S x {l') - 5,(2)5,(2') 

uc 

-5,(3)5,(3') + 5,(4)5,(4')] 
= N uc Q y Ky xx [t xx +t* xx ] , (C7) 

where 

t** = e- i<za [(a s , x +i6 SiX )(a* iX -i6* iX ) 
-(-a S;X + ib s , x ){-a* s x - i&*_ x ) 

- ib s,x)(a* s , x + ib *s, x ) 
+(-a s ,x - ^*,x)(-a*, x + iK,x)] 
= e-^[-4m s , x 6: iX +4ia: ;X 6 S;X ] 
= 8e-^S[a s , x 6; x ] , (C8) 

so that 

V^J = 16N uc Q y Ky xx co S (qa)S[a s , x bl x ] . (C9) 
Next, we consider terms proportional to -fQf,: 

V y (tf» y ) = Q,i^,^[5,(l)5,(l')- 5,(2)5,(2') 

-uc 

-5,(3)5,(30 + 5,(4)5,(4')] 
= N uc Q y Ky y [t yy + t; y ] , (CIO) 



where 



tyy = e ^"[("^ + b s,y)(-ia* s , y + K, y ) 
-{ia Si y - b St y)(-ia* s y - 6* y ) 



-(-ia S; , + b a>y )(ia* <y + &*,,) 
+ (-ia s> , - 6 s ,,)(m* y - 6* „)] 
= e-*«° [4ia a , B 6* w - 4»a* w 6 StS ] 
= 8e-^9[a: ;!/ 6 s ,,] , (Cll) 

so that 

V y (KV yy ) = 16N uc Q y Ky y cos(qa)Q[a* s J s J .(C12) 
Next, we consider terms proportional to K\ z : 
V v {Kl z ) = QyK\ z £(5,(1)5,(1') -5,(2)5,(2') 

UC 

-5,(3)5,(3') + 5,(4)5,(4')) 
= N uc Q y Ky zz [t zz +t* zz ] , (C13) 

where 

tzz = e ~ lqa [( a s,z+ib s , z )(a* z -ib* z ) 
-(a SiZ -ib s , z )(a* s z +ib* s z ) 
-(a StZ -ib StZ )(a* tZ + ib* z ) 
+(o SlZ + i6 s ,,)(a* j , - i&* J] 
= e-* qa [-4ia s , z bl z + 4ia* z 6 a ,»] 
= 8 e -^3[ as ,,6 s %] , (C14) 

so that 

V y {K\ z ) = 16N uc Q y Ky zz co S (qa)S[a s , z b* s J.(C15) 
Now we consider terms proportional to E^: 
V V {E«) = Q,^]T[5,(1)5,(1')- 5,(1)5,(1') 

UC 

+5,(2)5,(2')-5,(2)5,(2') 
+5,(3')5,(3)-5,(3')5,(3) 
+5,(4')5,(4)-5,(4')5,(4)] 
= N uc Q y Ey[t zy +t* zy } , (C16) 

where 

tz V = e~ lqa [(a s , x + tb s , x )Hal y + bl y ) 

-(ia s ,y + b s ,y)(a* s ,x- ib *s,x) 

+{~ a s,x + ib s ,x)(-ia* s , y - K,y) 

-(ia 8 , y - b Siy )(-a* tX - ib* s x ) 

+(a S;X - ib StX ){ia* s y +b* y ) 

-{-ia a ,y + b s ,y){a* >x + ib* tX ) 

+{-a SlX - ib^ x )(ia* y - b* y ) 

-{-ia s , y - b s , y )(-a* iX + ib* J] 
= 4e 9 [a S;X & SJ/ + b s , x a sy — &s,yb SiX — b Sl ,ci SjX ] 
= 8ie- lc,a Z{a s , x bl y + 6 s ,,ay , (C17) 

so that 

V V (E») = 16N uc Q y E y sm(qa)^s[a SyX b* s y + 6 SiX a*(pl8) 
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Now we consider terms proportional to 
V y {El) = Q v E»Y\S y {l)S t {l')- S Z (1)S V (1') 

UC 

-S y (2)S z (2') + S z (2)S y (2>) 
+S y (3')S z (3) - S z (3')S y (3) 
-S y (4')S z (4) + S z (4')S y (4)] 
= N uc Q y Ey[t xy +t* xy ] , (C19) 

where 

t xy = e- lqa [{ia s , y + b s , y ){al z ~ibl z ) 

-(a s ,z + ib s ,z)(-ia* s ,y + K, y ) 
-(ia s , y ~ b s , y )(a* tZ + ib* tZ ) 
H a s,z - ib s , z )(-ia* :y - b* y ) 
+(-ia S!y + b s , y )(a* iZ +ib* iZ ) 
-(a s , z -ib s , z )(ia* ty + bl y ) 
-(-ia s , y - b s , y )(a* z - ib* z ) 
+(a SiZ + ib s , z )(ia* !y - b* s y )] 

— Ap- i <i a \ n h* 4- h n* 

— [U s ^ y U s z T L>s,yU s z 

-0-s,zK,y - b s,* a *s,y] 

= Hie-^la^ + b^J , (C20) 

so that 

V y (E») = 16N uc Q y Ey S m(qa)Q[a s , y b* SiZ + b s , y a^21) 
Now we consider terms proportional to KJj. z : 
V y {Kl z ) = Q y Kl z Y,[S x {l)S z {l') + S z {l)S x {l') 

UC 

+S X {2)S Z (2') + S Z (2)S X {2') 
-S X (3)S Z (3') - S Z (3)S X (3') 
-5 X (4)S Z (4')-S,(4)5 X (4')] 
= N uc Q y Kl z [t xz +t* xz ] , (C22) 

where 

txz = e- iqa [{a s , x + ib s , x ){a* sz -ib* sz ) 

+{a s , z +ib s , z )(a* s x - ib* s x ) 

+(-a SiX + ib s , x ){a* z + ib* z ) 

+{a s , z -ib StZ ){-a s x -ib* s x ) 
-{a s , x -ib SiX )(a* tZ + ib* z ) 

-(a S;2 -ib s , z )(a* s x + ib* s x ) 

-{-a SlX - ib s , x )(a* s z - ib* z ) 

-(a S;2 + ib s , z )(-a* sy + ib* sy )] 
= 4:ie~ tqa [-a SyX b* z + b StX a* z - a s , z b* x + b StZ a* tX ] 
= 8e-* qa Z{a s , x bl z + a s , z b* s J , (C23) 

so that 

V v (K% t ) = l6N uc Q y Ky xz cos(qa) 

x%[a s , x bl z + a s , z b* s J . (C24) 



3. Spine cross-tie interactions 

Here we analyze the spin-cross tie interactions. 
We consider the terms proportional to L xx : 

V y (LV x ) = Q y Ll x [S x (5) + S x (6)] 

x [S x (l) -S x (2) -S x (3) + S x (4)] 



— N uc Q y L xx [u xx + u x 



(C25) 



where 

u xx = b cx [(a* x ~ib* x )e-' iqa,A 



-(-<* - K ^ iqa/i - K* + K*Y qa,i 

+(-<* + K*y qa/i 

-«* - K*Y qa,i + (-<* K*Y qa/i 

+«* + rt*)^ iqa,i - (-<* + ib* )e~ iqa/i ] 



= 4b cx a*Je' lqa /' i - e lqa/i ] 



(C26) 



= —&ism{qa/A)b cx a* x 
so that 

Vy{V xx ) = 16N uc QyLy xx sm(qa/4)S[b cx a* SiX ] (C27) 
We consider the terms proportional to L v yy : 

Vy(Pyy) = QyL V yy [S y (5) + S y (6)} 

[S y (l) - S y (2) - S y (3) + S y (4)] 



— N U cQyL yy [ u yy 



(C28) 



where 

U yy — Q-cy[(~ ^ a s,y b s .y) e ^ ^ 



(-K, y - bt y )e-^ 
-«, v + K, y )e iqa ' 4 + (ial y - 6:, y ) e ^/ 4 
-(-<„ + bl y y qa / 4 + (-ia* Sty - bl y )e^ 
+«,y + bt y )e-* qa/i «„ - KJe-^} 
= Aa cy bl y [e- tqa ^ - e^/ 4 ] 

= -8ism(qa/4)a cy b* y (C29) 
so that 

V v {Ly y ) = 16N uc Q y Ly yS m(qa/4)Z[a cy b* s J .(C30) 

We consider the terms proportional to L v zz : 

V y {Ll z ) = QyLl z [S z {h) + S z {Q)] 

x[S z (l) - S z (2) - S z (3) + S z (4)] 
= N uc Q y Ly zz [u zz +u* zz ] , (C31) 

where 

u zz = a cz [(a* StZ - l KJe-^ i « z + ^ Je"^ 4 
H< z + KzV qa/i + Kz-*bizV va/i 

+«, z KzY qa,i ~ Kz + *KzV qa/i 
-K* + Kz)^ tqa/i + Kz Kz)^ tqa/i ] 



4a cz b* )Z i 



_ e -i<?a/4 _ e iqa/4 



= -8ia cz b* z cos(aq/4) , 



(C32) 
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so that 

V V (LV X ) = 16N uc Q y Ly zz co S (qa/4)Z[a cz b* s J .(C33) 
We consider the terms proportional to L v xy : 

V y (L% y ) = Q v L%([S x {h)-S x (<S)] 

X[Sy(l) + Sy(2) - Sy(3) - Sy(4)] 

+[S V (5) 5,(6)] 

x[S x (l) + S x (2) - S x (3) - S x (4)} 
= N uc Q v Ll y [u xy + u* xy ] , (C34) 



where 



{-i< v +b*y- iqa/A +{- 



10 



-iqa/A 



^< v + K. v y qa/i -^^-b:„u 



s,y "s,y 
iqa/A 



s,yJ 



+(-< v + K, y )e iga/4 + (-<„ - K, y )e iga/4 



cy 



— n* p~ ic l a l^ — n * p'W 4 

s,y s,y 



~K x -iKJe- i9a/4 

+(-< x - K^ iqa,i - «* + K*y qa,i 
-(-<, x +K*y qa/i 

-«* + iblJe-W* (_ * a + ^Je"^ 4 

= -8i[6 CK a* !y +a C j / 6* x ]cos(ga/4) , (C35) 
so that 

V y (Ll y ) = 16N uc Q v Ly y cos(qa/4) 

x%[b cx al y + a cy bl x ] . (C36) 

We consider the terms proportional to L xz : 

V V (L* X ) = Q y Ll z ([S x (5) + S x (6)) 

x[S z (l) + S z (2) - S z (3) - S z (4)} 
+ [S Z (5) + S z (6)} 

x[S x (l) + S x {2) - S x (3) - S X (4)]J 
= N uc Q y Ll z [u xz + u* xz ] , (C37) 



where 



Kz < z )e- iqa/4 + Kz + Kz)e- iqa/4 



-Kz + KzY qa ' 4 Kz KzY qa/4 



-Kz iK,z)e iqa ' 4 Kz + KzY qa/4 

+Kz + Kz)e- iqa/i + Kz - Kz)e- iqa/4 
+a cz [« x - iKJe-^ 

+(-<* - K*) e ~ iqa/4 

-«* + <*Y qa/A (-<* + < x Y qa/4 

+ K X ibl x )e iqa ' A + (-< x - iK tX )eW 

-«* + K x )e- iqa/4 (-«:,, + iK, x )e-^ a / 4 

= 4b cx a* s Je- lqa/4 - e - lqa/i ] 

+4a cz b*J-ie~ tqa/4 -ie tqa/A ] 
= —8ib cx al z sm(qa/4) — 8ia cz b* x cos(qa/4) ,(C38) 

so that 

V v {Ll z ) = 16N uc Q y Ly xz ( S m(qa/^[b cx a* StZ } 



+ cos(qa/4)Z[a cz b* s jj . (C39) 
We consider the terms proportional to L v yz : 

V V (LV X ) = g,L»,|[S„(5)-S,(6)] 

x[S z (l) - S z (2) - S z (3) + S z (4)} 
+ [S Z (5) - S z (6)} 

X[Sy(l)-Sy(2)-Sy(3)+Sy(4)]j 

= N uc QyLy yz [uy Z + u* vz ] , (C40) 



where 



' 1/2 



•■cy 



-( a : !Z+J 6: ! je^/ 4 + (a: iZ - J 6: i je^/ 4 

+«, - ^;je^/ 4 - « 2 + i& s %)e^/ 4 
-Kz + Kz>~ iqa/4 + Kz - Kz>~ iqa/4 

+a c , z [(-ialy + ft* Je"* / 4 

-(-Kv-Ky)e- iqa/4 

~(Kv + KyY qa/i + (Kv ~ KyY^ 4 

-(-*<» + & ;> iW4 + (-<„ - &:, y )^ a/4 



— 4ia cy b sz 



-4a cz b* y 



e -iga/4 giga/4 



g -iga/4 _ iqa/i 



-8ia cy b* s z cos(aq/4) - 8ia cz b* y sin(ga/4) ,(C41) 
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so that 

Vy(Ly Z ) = 16N uc Q y Ly yz (cos(qa/4)Z[a C yb* s J 

+ sm(qa/4)S[a cz b*J) . (C42) 

We consider the terms proportional to F%: 

Vy(Fy) = Q y F%([S y (5)-S y (6)} 

x [S z (l) -S z (2) -5,(3) + S Z (4)\ 
~[S Z (5) - S z (6)\ 

x[S y (l)-S y (2)-S y (3) + S y (4)]j 
= N uc Q y Fy[u x + u* x ] , (C43) 



where 



where 



•■cy 



+«, - K zY qa/4 - K* + KzY qa/4 

+ K z)e-' qa/i + (a* - ib* )e-^ 



+a c 



-iqa/4 



V 1 s,yJ 
iqa/4 



+ (- ia *s,y - b ty)e 

+{ia* Sty + b* V'W4_ 



iqa/4 



) e -iga/4 + (m 



s,y ®s,y)^ 



-iqa/A 



Aia C yb s z 



+4a cz b* s>y 



-iqa/A , & iqa/A 



e iqa/4 e -i(ja/4 



-$ia C yb* s z cos(qa/4) + 8m cz 6* ;J/ sin(ga/4) ,(C44) 



so that 



- l6N uc Q y Fv[cos(qa/4:)Z[a cy b* s J 

+ S m(qa/4)<Z[a* cz b s , y } ) . (C45) 



We consider the terms proportional to F^jj: 

Vy(Fy) = Q y Fy(-[S x (5) + S x (6)} 

x [S z (l) + S Z (2) -S z (3) -S z (4)\ 
+ [S Z (5) + 5,(6)] 

x[5 x (l) + 5 x (2)-5 x (3)-5 x (4)fj 
= N uc Q y Ly yz [u y + u* y ] , (C46) 



+«, + KzY qa/i + «z - KzY qa/i 
H<,z-KzY qa/i + {<, z + iK,zY qa/i 

-Kz + Kz)^ tqa/i - Kz Kz)^ tqa/i 



-iqa/4 



+(-<*- K*) e ~ iqa/4 

-«* + K*Y qa,i ("<x + ^:,x)e <9o/4 

«* - j^ a / 4 + (-< x - J 6:,j e ^/ 4 



+4a cz 6! 



-je^ a/4 - ie l9Q/4 
= 8ib cx a* z sin(ga/4) — 8ia cz b* x cos(ag/4) , (C47) 

so that 



K,(^) = lW uc Q y FV y Un{qa/4)%b* cx a StZ 



■cos(aq/4)$s[a cz b* x 



(C48) 



Finally, we consider the terms proportional to F^: 
Vy(Fy) = Q y Fy([S x (5)-S x (6)} 

X[Sy(l) + Sy(2) - Sy(3) - Sy(4)} 

-[5,(5) -5,(6)] 

x[5 x (l) + 5 x (2)-5 x (3)-5 x (4)]^ 
= N uc Q y Fy[u z +u* z ] , (C49) 



where 

u, = b 



H-Kv + b*Y qa/4 + H< 



iqa/4 



■ KyY 

-«v + Kv>~ iqa/i - (K « - K Je"^ 4 



+a. 



c,y 



y 1 s-v> \ s,y 

qa/4 



(a* s ,x ~ ib *s,xY 

-(-<* - K, x )z- iqa '* 

+«* + iK, x Y qa/i + (-<* + ^:,x)e <9o/4 
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-iqa/4 



so that 



+Ua cy b*^ x 



Vy(FV) = 16iV„ c Q H F|cos(ga/4) 



x3[6 a a: !V + a;6 8iI ]. (C51) 



-8?6 CT a* cos(ga/4) + 8ia cy b* s x cos(ga/4) , (C50) 
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